clear
est clear

	cap mat define migrationfigure=J(105,10,.)
		local row=1
		

	use "$input/full_geoid_pair_panel.dta"

			
	replace migration=migration/(tot_pop/100000)

	//Make sure only include states that have a storm in at least one county at least once in our sample years
	bysort state_fips: egen max=max(storm)
		drop if max!=1
		
	//Create future pooled migration - total migrants in the subsequent five years
	sort geoid_pair move_year
	g post_storm=1 if storm==1|storm[_n-1]==1&geoid_pair==geoid_pair[_n-1]|storm[_n-2]==1&geoid_pair==geoid_pair[_n-2]|storm[_n-3]==1&geoid_pair==geoid_pair[_n-3]|storm[_n-4]==1&geoid_pair==geoid_pair[_n-4]|storm[_n-5]==1&geoid_pair==geoid_pair[_n-5]
	mvencode post_storm, mv(0)
	
	//Collapse to total migration from a county (not looking at paired county to county migration, just all migration out of a county)
	collapse (sum) migration (mean) share_approved totalapprovedihpamount,by(from_geoid move_year storm storm_fema post_storm state_fips)
	
	//Create IHS transformations
	ihstrans(migration)
	
	g total_storms=0
	//Create lags of storms (was there a storm in t-1 years)
	sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g storm_L`i'=storm[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
				replace total_storms=total_storms+storm[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
			
	label var post_storm "1(0-5 yrs post storm)"
	label var total_storms "All storms in t-5 years"
	label var storm "1(Storm year)"
	foreach var in 1 2 3 4 5{
		label var storm_L`var' "1(Storm year), t-`var'"
	}
	
	//Create leads of storms (was there a storm in t+1 years)
	sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g storm_E`i'=storm[_n+`i'] if from_geoid==from_geoid[_n+`i'] 
				 
		}	
			
	
	foreach var in 1 2 3 4 5{
		label var storm_E`var' "1(Storm year), t+`var'"
	}
	
	
	
	
	replace totalap=totalap/1000000
	
		
		reghdfe ihs_migration i.storm, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			
			
			
			eststo t1c1
			cap mat def migrationfigure[`row',1]=_b[1.storm]
			cap mat def migrationfigure[`row',2]=_se[1.storm]
			cap mat def migrationfigure[`row',3]="no lags, out"
			g sample=e(sample)
			local row=`row'+1
			
				
		reghdfe ihs_migration i.storm i.storm_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t1c2
			cap mat def migrationfigure[`row',1]=_b[1.storm]
			cap mat def migrationfigure[`row',2]=_se[1.storm]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L1]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L1]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L2]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L2]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L3]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L3]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L4]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L4]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L5]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L5]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			
			
			lincom 1.storm+1.storm_L1+1.storm_L2+1.storm_L3+1.storm_L4+1.storm_L5
			
				local beta11=trim("`: display %10.3f r(estimate)'")
				local se11=trim("`: display %10.3f r(se)'")	
				cap mat def migrationfigure[`row',1]=`beta11'
				cap mat def migrationfigure[`row',2]=`se11'
				cap mat def migrationfigure[`row',3]="lags, out"
				local row=`row'+1
				
		reghdfe ihs_migration i.post_storm, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			*eststo t1c3
				
		
		reghdfe ihs_migration i.storm i.storm_E*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c1
			
		reghdfe ihs_migration i.storm i.storm_E* i.storm_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c2
		
		reghdfe ihs_migration i.storm total_storms, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t3c1
		
		reghdfe ihs_migration i.storm##c.total_storms, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t3c2
			
		keep from_geoid sample storm
		drop if sample!=1
		bysort from_geoid: egen total_storms=total(storm)
		drop storm
		duplicates drop from_geoid, force
		rename from_geoid fips
		export delimited "$figdat/migration_sample.csv", replace
	
	clear
	
	use "$input/full_geoid_pair_panel.dta"
	
	bysort state_fips: egen max=max(storm)
		drop if max!=1
		
	sort geoid_pair move_year
	g post_storm=1 if storm==1|storm[_n-1]==1&geoid_pair==geoid_pair[_n-1]|storm[_n-2]==1&geoid_pair==geoid_pair[_n-2]|storm[_n-3]==1&geoid_pair==geoid_pair[_n-3]|storm[_n-4]==1&geoid_pair==geoid_pair[_n-4]|storm[_n-5]==1&geoid_pair==geoid_pair[_n-5]
	mvencode post_storm, mv(0)	
		
	replace migration=migration/(tot_pop/100000)
	
			
	preserve
		collapse (sum) migration ,by(to_geoid move_year )

			rename migration in_migration
			rename to_geoid from_geoid
			save "$input/temp1.dta", replace
	restore
	
	
		collapse (sum) migration ,by(from_geoid move_year storm post_storm state_fips)
	
		merge 1:1 from_geoid move_year using "$input/temp1.dta"
		
		g net_migration=in_migration-migration
		
		sort from_geoid move_year
		//Create lags of storms (was there a storm in t-1 years)
		g total_storms=0
		foreach i in 1 2 3 4 5{
					g storm_L`i'=storm[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
					replace total_storms=total_storms+storm[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
			
		label var post_storm "1(0-5 yrs post storm)"
		label var total_storms "All storms in t-5 years"
				
		
		label var storm "1(Storm year)"
		foreach var in 1 2 3 4 5{
			label var storm_L`var' "1(Storm year), t-`var'"
		}
		
		//Create leads of storms (was there a storm in t+1 years)
		sort from_geoid move_year
			foreach i in 1 2 3 4 5{
					g storm_E`i'=storm[_n+`i'] if from_geoid==from_geoid[_n+`i'] 
					 
			}	
				
		
		foreach var in 1 2 3 4 5{
			label var storm_E`var' "1(Storm year), t+`var'"
		}
			
		ihstrans(net_migration)	
		
		
		
		

		reghdfe ihs_net_migration i.storm, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			
			
			eststo t1c4
			cap mat def migrationfigure[`row',1]=_b[1.storm]
			cap mat def migrationfigure[`row',2]=_se[1.storm]
			cap mat def migrationfigure[`row',3]="no lags, net"
			local row=`row'+1
		
				
		reghdfe ihs_net_migration i.storm i.storm_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t1c5
			cap mat def migrationfigure[`row',1]=_b[1.storm]
			cap mat def migrationfigure[`row',2]=_se[1.storm]
			cap mat def migrationfigure[`row',3]="lags, net"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L1]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L1]
			cap mat def migrationfigure[`row',3]="lags, net"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L2]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L2]
			cap mat def migrationfigure[`row',3]="lags, net"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L3]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L3]
			cap mat def migrationfigure[`row',3]="lags, net"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L4]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L4]
			cap mat def migrationfigure[`row',3]="lags, net"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.storm_L5]
			cap mat def migrationfigure[`row',2]=_se[1.storm_L5]
			cap mat def migrationfigure[`row',3]="lags, net"
			local row=`row'+1
			lincom 1.storm+1.storm_L1+1.storm_L2+1.storm_L3+1.storm_L4+1.storm_L5
			
				
				
		reghdfe ihs_net_migration i.post_storm, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			*eststo t1c6
				
		reghdfe ihs_net_migration i.storm i.storm_E*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c5	
		
		reghdfe ihs_net_migration i.storm i.storm_E* i.storm_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c6	
			
		reghdfe ihs_net_migration i.storm total_storms, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t3c4
		
		reghdfe ihs_net_migration i.storm##c.total_storms, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t3c5
			
		
	graph set window fontface Arial
	preserve
		clear
			
		svmat2 migrationfigure

		drop if migrationfigure1==.
	
		g n=_n
			replace n=1.8 if n==2
			replace n=2.4 if n==3 
			replace n=3 if n==4
			replace n=3.6 if n==5
			replace n=4.2 if n==6
			replace n=4.8 if n==7
			replace n=5.6 if n==8
		
		g HI95_mean=(1.96*migrationfigure2)+migrationfigure1
		g LI95_mean=migrationfigure1-(1.96*migrationfigure2)

		g HI99_mean=(2.06*migrationfigure2)+migrationfigure1
		g LI99_mean=migrationfigure1-(2.06*migrationfigure2)
		
		set scheme plotplain
		
		rename migrationfigure1 mean_beta
		
		replace mean_beta=mean_beta*100 
		replace HI99_mean=HI99_mean*100 
		replace LI99_mean=LI99_mean*100
		
		replace n=n-1
	
		twoway (bar mean_beta n  if n<8, barwidth(0.5) color("0 45 114%50")) ///
					(rspike HI99_mean LI99_mean n if n<8, lc(gs8)  ///
					yti("Percent change in out migration",size(vlarge)) ylabel(, nogrid labsize(large)) xlabel(0 "Storm year" 0.8 "Storm year" 1.4 "Storm year, L1" 2 "Storm year, L2" 2.6 "Storm year, L3" 3.2 "Storm year, L4" 3.8 "Storm year, L5" 4.6 "Sum of lags",nogrid labsize(large)) ysize(1.5) xsize(5) xtitle("") legend(off) yline(0) xline(0.4) xline(4.2))
					
					graph export "$figures/fig_1c.png", as(png) replace
					
					
	restore
		
	preserve
		clear
			
		svmat2 migrationfigure

		drop if migrationfigure1==.
		
		drop in 1/8
	
		g n=_n
			replace n=1.8 if n==2
			replace n=2.4 if n==3 
			replace n=3 if n==4
			replace n=3.6 if n==5
			replace n=4.2 if n==6
			replace n=4.8 if n==7
			replace n=5.6 if n==8
		
		g HI95_mean=(1.96*migrationfigure2)+migrationfigure1
		g LI95_mean=migrationfigure1-(1.96*migrationfigure2)

		g HI99_mean=(2.06*migrationfigure2)+migrationfigure1
		g LI99_mean=migrationfigure1-(2.06*migrationfigure2)
		
		set scheme plotplain
		
		rename migrationfigure1 mean_beta
		
		replace mean_beta=mean_beta*100 
		replace HI99_mean=HI99_mean*100 
		replace LI99_mean=LI99_mean*100
		
		replace n=n-1
	
		twoway 	(bar mean_beta n  if n<8, barwidth(0.5) color("252 89 16%50")) ///
					(rspike HI99_mean LI99_mean n if n<8, lc(gs8)  ///
					yti("Percent change in net migration",size(vlarge)) ylabel(, nogrid labsize(large)) xlabel(0 "Storm year" 0.8 "Storm year" 1.4 "Storm year, L1" 2 "Storm year, L2" 2.6 "Storm year, L3" 3.2 "Storm year, L4" 3.8 "Storm year, L5" 4.6 "Sum of lags",nogrid labsize(large)) ysize(1.5)  xsize(5) xtitle("") legend(off) yline(0) xline(0.4) xline(4.2))
					
					graph export "$figures/fig_1d.png", as(png) replace
					
					
	restore
	
			
	